Cost of facilitation: invasive grasses limit recruitment of benefactor shrubs.

ecoblender

alex.filazzola

Abstract

Shrubs facilitate the abundance and productive of annual plants in desert ecosystems.However, these shrub microhabitats favour plant species with competitive life histories. Ephedra californica is a dominant shrub in the San Joaquin desert that has been identified as a facilitator. Here, we explore the factors that limit the Ephedra californica recruitment into the San Joaquin desert including substrate, water availability, and herbivory. We also explore the role of the invasive grass Bromus madritensis on limiting establishment of E. californica. Vegetation surveys were conducted in the field during 2013 and collected seed was then used to conduct two greenhouse trials. The first explore germination and establishment techniques for E. california and the optimal substrate. The second examined E. californica establishment in present of the invasive B. madritensis responding to different water levels and herbivory. These results can have implications for land managers in the San Joaquin Valley to maintain native shrub biodiversity in the region.

---

## load libraries
library(tidyverse)
library(OIsurv)
library(lsmeans)
library(ggthemes)

## load data
substrate <-read.csv("data/Ephedra.substrate.csv")
recruit <-read.csv("data/Ephedra.recruitment.csv")
landscape <-read.csv("data/ephedra.landscape.csv")

##load functions
## inverse hyperbolic sine transformation 
ihs <- function(x) {
    y <- log(x + sqrt(x ^ 2 + 1))
    return(y)
}
se <- function(x, ...) sqrt(var(na.omit(x))/length(na.omit(x)))
source("functions.r")

Ephedra characteristics at the landscape

avg.size <- landscape %>% summarize(H=mean(H), D1=mean(D1), Area=mean(Area), density=mean(Shrub.density), rdm=mean(RDM.2013))
avg.size

##          H       D1     Area     density      rdm
## 1 1.332925 3.446412 9.321902 0.004318859 7.042223

Ephedra recruitment at landscape

## compare ephedra size and density across landscape
landscape[,"area.group"] <- as.numeric(cut(landscape$Log.area,20))
hist(landscape$area.group)

mean.area <- landscape %>% group_by(area.group) %>%  summarize(rdm=mean(RDM.2013), area=mean(Log.area), density=mean(Shrub.density))

ggplot(mean.area) + geom_point(aes(x=density, y=area), size=3)+ylab("log (area of shrub)")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab(expression("shrub density (m"^"2"*")")) +  stat_smooth(method="lm", formula= y~x,aes(x=density, y=area))

m1 <- lm(area~density, data=mean.area)
summary(m1)

## 
## Call:
## lm(formula = area ~ density, data = mean.area)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.82510 -0.26068  0.01197  0.22216  1.16732 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     9.3959     0.9852   9.537 1.84e-08 ***
## density     -1706.0323   223.6888  -7.627 4.81e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5569 on 18 degrees of freedom
## Multiple R-squared:  0.7637, Adjusted R-squared:  0.7506 
## F-statistic: 58.17 on 1 and 18 DF,  p-value: 4.806e-07

ggplot(mean.area) + geom_point(aes(x=density, y=rdm), size=mean.area$area+1)+ylab("residual dry matter")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab(expression("shrub density (m"^"2"*")")) +  stat_smooth(method="lm", formula= y~exp(x),aes(x=density, y=rdm))

m2 <- lm(log(rdm)~density, data=mean.area)
summary(m2)

## 
## Call:
## lm(formula = log(rdm) ~ density, data = mean.area)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32696 -0.09127 -0.02044  0.05849  0.46804 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.1851     0.3294   9.670 1.49e-08 ***
## density     -279.3744    74.7892  -3.735  0.00151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1862 on 18 degrees of freedom
## Multiple R-squared:  0.4367, Adjusted R-squared:  0.4054 
## F-statistic: 13.95 on 1 and 18 DF,  p-value: 0.001514

Ephedra optimal substrate

substrate[,"Sand"] <- as.factor(substrate$Sand)
substrate[,"census"] <- as.factor(substrate$census)

m1 <- aov(height ~ Micro * Sand,  data=subset(substrate, census==10))
summary(m2) ## nothing significant

## 
## Call:
## lm(formula = log(rdm) ~ density, data = mean.area)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32696 -0.09127 -0.02044  0.05849  0.46804 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.1851     0.3294   9.670 1.49e-08 ***
## density     -279.3744    74.7892  -3.735  0.00151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1862 on 18 degrees of freedom
## Multiple R-squared:  0.4367, Adjusted R-squared:  0.4054 
## F-statistic: 13.95 on 1 and 18 DF,  p-value: 0.001514

## compare shade on survival of ephedra
my.surv <- Surv(as.numeric(substrate$census), substrate$survival)
fit1 <- survfit(my.surv~Micro, data=substrate)
summary(fit1)

## Call: survfit(formula = my.surv ~ Micro, data = substrate)
## 
##                 Micro=Open 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    960      11    0.989 0.00343        0.982        0.995
##     2    880      29    0.956 0.00681        0.943        0.969
##     3    800      41    0.907 0.00987        0.888        0.927
##     4    720      42    0.854 0.01221        0.830        0.878
##     5    640      67    0.765 0.01504        0.736        0.795
##     6    481      31    0.715 0.01647        0.684        0.748
##     7    400      62    0.604 0.01901        0.568        0.643
##     8    241      31    0.527 0.02108        0.487        0.570
##     9    160      31    0.425 0.02366        0.381        0.474
##    10     80      31    0.260 0.02730        0.212        0.320
## 
##                 Micro=Shade 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    960      33    0.966 0.00588        0.954        0.977
##     2    880      43    0.918 0.00897        0.901        0.936
##     3    800      42    0.870 0.01117        0.849        0.892
##     4    720      40    0.822 0.01290        0.797        0.848
##     5    640      34    0.778 0.01422        0.751        0.807
##     6    480      17    0.751 0.01521        0.721        0.781
##     7    400      31    0.692 0.01725        0.659        0.727
##     8    240      16    0.646 0.01958        0.609        0.686
##     9    160      15    0.586 0.02317        0.542        0.633
##    10     80      13    0.491 0.03099        0.433        0.555

par(mar=c(4.5,4.5,.5,.5))
plot(fit1, col="white", ylim=c(0.2,1), xlim=c(0.9,10.1), ylab="Estimated Survival Function", xlab="Census", cex.lab=1.5, cex.axis=1.3)
lines(1:10,summary(fit1)[[10]][1:10], col="#E69F00", lty=2) ## upper
lines(1:10,summary(fit1)[[6]][1:10], col="#E69F00", lwd=2) ## value sun
lines(1:10,summary(fit1)[[11]][1:10], col="#E69F00", lty=2) ## lower
lines(1:10,summary(fit1)[[10]][11:20], col="#56B4E9", lty=2) ## upper
lines(1:10,summary(fit1)[[6]][11:20], col="#56B4E9", lwd=2) ## value shade
lines(1:10,summary(fit1)[[11]][11:20], col="#56B4E9", lty=2) ## lower
legend(1.8, 0.35, c("Sun","Shade"), lty=1, lwd=3, col=c("#E69F00","#56B4E9"), cex=1.5)

coxph.fit <- coxph(my.surv ~ Micro, method="breslow", data=substrate)
coxph.fit 

## Call:
## coxph(formula = my.surv ~ Micro, data = substrate, method = "breslow")
## 
##               coef exp(coef) se(coef)     z       p
## MicroShade -0.2802    0.7557   0.0786 -3.56 0.00037
## 
## Likelihood ratio test=12.8  on 1 df, p=0.000342
## n= 1920, number of events= 660

## compare sand on survival of ephedra
my.surv <- Surv(as.numeric(substrate$census), substrate$survival)
fit2 <- survfit(my.surv~Sand, data=substrate)
summary(fit2)

## Call: survfit(formula = my.surv ~ Sand, data = substrate)
## 
##                 Sand=0 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       7    0.982 0.00683        0.968        0.995
##     2    352      11    0.951 0.01125        0.929        0.973
##     3    320      13    0.912 0.01506        0.883        0.942
##     4    288      12    0.874 0.01799        0.840        0.910
##     5    256      14    0.827 0.02106        0.786        0.869
##     6    192       6    0.801 0.02289        0.757        0.847
##     7    160      12    0.741 0.02695        0.690        0.795
##     8     96       6    0.694 0.03120        0.636        0.758
##     9     64       6    0.629 0.03794        0.559        0.708
##    10     32       6    0.511 0.05325        0.417        0.627
## 
##                 Sand=25 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384      10    0.974 0.00813        0.958        0.990
##     2    352      16    0.930 0.01331        0.904        0.956
##     3    320      20    0.872 0.01772        0.838        0.907
##     4    288      18    0.817 0.02075        0.777        0.859
##     5    256      15    0.769 0.02292        0.726        0.815
##     6    192       7    0.741 0.02441        0.695        0.791
##     7    160      14    0.676 0.02776        0.624        0.733
##     8     96       7    0.627 0.03137        0.568        0.692
##     9     64       7    0.558 0.03714        0.490        0.636
##    10     32       7    0.436 0.05007        0.348        0.546
## 
##                 Sand=50 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384      11    0.971 0.00851        0.955        0.988
##     2    352      16    0.927 0.01350        0.901        0.954
##     3    320      19    0.872 0.01765        0.838        0.907
##     4    288      20    0.812 0.02098        0.771        0.854
##     5    256      18    0.755 0.02343        0.710        0.802
##     6    192       9    0.719 0.02512        0.672        0.770
##     7    160      17    0.643 0.02848        0.589        0.701
##     8     96       8    0.589 0.03178        0.530        0.655
##     9     64       8    0.516 0.03697        0.448        0.593
##    10     32       8    0.387 0.04823        0.303        0.494
## 
##                 Sand=75 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       9    0.977 0.00772        0.962        0.992
##     2    352      18    0.927 0.01361        0.900        0.954
##     3    320      16    0.880 0.01716        0.847        0.915
##     4    288      17    0.828 0.02025        0.790        0.869
##     5    256      28    0.738 0.02422        0.692        0.787
##     6    192      13    0.688 0.02624        0.638        0.741
##     7    160      26    0.576 0.02976        0.521        0.637
##     8     96      13    0.498 0.03266        0.438        0.566
##     9     64      13    0.397 0.03612        0.332        0.474
##    10     32      13    0.236 0.04058        0.168        0.330
## 
##                 Sand=100 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    384       7    0.982 0.00683        0.968        0.995
##     2    352      11    0.951 0.01125        0.929        0.973
##     3    320      15    0.907 0.01554        0.877        0.937
##     4    288      15    0.859 0.01891        0.823        0.897
##     5    256      26    0.772 0.02349        0.727        0.819
##     6    193      13    0.720 0.02596        0.671        0.773
##     7    160      24    0.612 0.03000        0.556        0.674
##     8     97      13    0.530 0.03351        0.468        0.600
##     9     64      12    0.431 0.03755        0.363        0.511
##    10     32      10    0.296 0.04372        0.222        0.395

test <- substrate %>% group_by(census,Sand) %>%  summarize(count=sum(survival),avg=mean(survival))
test <- data.frame(test)

par(mar=c(4.5,4.5,.5,.5))
plot(fit2, col="white", ylim=c(0.2,1), xlim=c(0.9,10.1), ylab="Estimated Survival Function", xlab="Census", cex.lab=1.5, cex.axis=1.3)
lines(1:10,summary(fit2)[[10]][1:10], col="#E69F00", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][1:10], col="#E69F00", lwd=2) ## value sand 0
lines(1:10,summary(fit2)[[11]][1:10], col="#E69F00", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][11:20], col="#56B4E9", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][11:20], col="#56B4E9", lwd=2) ## value sand 25
lines(1:10,summary(fit2)[[11]][11:20], col="#56B4E9", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][21:30], col="#009E73", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][21:30], col="#009E73", lwd=2) ## value sand 50
lines(1:10,summary(fit2)[[11]][21:30], col="#009E73", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][31:40], col="#CC79A7", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][31:40], col="#CC79A7", lwd=2) ## value sand 75
lines(1:10,summary(fit2)[[11]][31:40], col="#CC79A7", lty=2) ## lower
lines(1:10,summary(fit2)[[10]][41:50], col="#D55E00", lty=2) ## upper
lines(1:10,summary(fit2)[[6]][41:50], col="#D55E00", lwd=2) ## value sand 100
lines(1:10,summary(fit2)[[11]][41:50], col="#D55E00", lty=2) ## lower
legend(1.7, 0.44, c("0%","25%","50%","75%","100%"), lty=1, lwd=3, col=c("#E69F00","#56B4E9","#009E73","#CC79A7","#D55E00"), cex=1.2, title="Sand")

coxph.fit <- coxph(my.surv ~ Sand, method="breslow", data=substrate)
coxph.fit 

## Call:
## coxph(formula = my.surv ~ Sand, data = substrate, method = "breslow")
## 
##          coef exp(coef) se(coef)    z       p
## Sand25  0.263     1.301    0.138 1.91  0.0563
## Sand50  0.365     1.441    0.135 2.71  0.0068
## Sand75  0.579     1.785    0.130 4.47 7.7e-06
## Sand100 0.450     1.568    0.133 3.39  0.0007
## 
## Likelihood ratio test=23.3  on 4 df, p=0.000109
## n= 1920, number of events= 660

## raw number of germinants end trials

end <- subset(substrate, census==10)

## compare microsite

end.means <- end %>% group_by(Micro) %>% summarize(eph=mean(survival), eph.se=se(survival), Height=mean(height, na.rm=T), Height.se=se(height)) %>% gather(measure, value, 2:5) %>% separate(measure, c("Eph",".se")) ## calculate metrics and put long format

## Warning: Too few values at 4 locations: 1, 2, 5, 6

end.means[is.na(end.means$.se),3] <- "mean"
end.means <- end.means %>% spread(.se, value)
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"mean"]-end.means[,"se"]
end.means[,"upper"] <- end.means[,"mean"]+end.means[,"se"]

survival <- subset(end.means, Eph=="eph")

plot1 <- ggplot(survival, aes(x=Micro, y=mean)) + geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Microsite") 

height <- subset(end.means, Eph=="Height")

plot2 <- ggplot(height, aes(x=Micro, y=mean)) + geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average height of Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Microsite") 

require(gridExtra)

## Loading required package: gridExtra

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

grid.arrange(plot1, plot2, ncol=2)

m1 <- glm(survival~ Sand* Micro, data=end, family="binomial")
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: survival
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                         159     188.21            
## Sand        4   4.7398       155     183.47 0.315056   
## Micro       1  10.7288       154     172.75 0.001055 **
## Sand:Micro  4   8.6588       150     164.09 0.070217 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lsmeans(m1, pairwise~Sand)

## NOTE: Results may be misleading due to involvement in interactions

## $lsmeans
##  Sand     lsmean          SE df   asymp.LCL    asymp.UCL
##  0    -1.4663371   0.4529108 NA   -2.354026  -0.57864819
##  25   -1.3671838   0.4643107 NA   -2.277216  -0.45715153
##  50   -8.7830342 494.5226042 NA -978.029528 960.46345952
##  75   -0.4236489   0.3831780 NA   -1.174664   0.32736619
##  100  -0.8047190   0.3872983 NA   -1.563810  -0.04562817
## 
## Results are averaged over the levels of: Micro 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast    estimate          SE df z.ratio p.value
##  0 - 25   -0.09915331   0.6486236 NA  -0.153  0.9999
##  0 - 50    7.31669717 494.5228116 NA   0.015  1.0000
##  0 - 75   -1.04268814   0.5932568 NA  -1.758  0.3987
##  0 - 100  -0.66161811   0.5959263 NA  -1.110  0.8012
##  25 - 50   7.41585049 494.5228222 NA   0.015  1.0000
##  25 - 75  -0.94353482   0.6020048 NA  -1.567  0.5185
##  25 - 100 -0.56246480   0.6046358 NA  -0.930  0.8852
##  50 - 75  -8.35938531 494.5227526 NA  -0.017  1.0000
##  50 - 100 -7.97831529 494.5227559 NA  -0.016  1.0000
##  75 - 100  0.38107003   0.5448168 NA   0.699  0.9567
## 
## Results are averaged over the levels of: Micro 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 5 estimates

## survival
end.means <- end %>% group_by(Sand) %>% summarize(eph=mean(survival), eph.se=se(survival))
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"eph"]-end.means[,"eph.se"]
end.means[,"upper"] <- end.means[,"eph"]+end.means[,"eph.se"]

ggplot(end.means, aes(x=Sand, y=eph)) + geom_bar(stat="identity", fill="#56B4E9") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") 

## height

end.means <- end %>% group_by(Sand) %>% summarize(eph=mean(height, na.rm=T), eph.se=se(height, na.rm=T))
end.means <- data.frame(end.means)

end.means[,"lower"] <- end.means[,"eph"]-end.means[,"eph.se"]
end.means[,"upper"] <- end.means[,"eph"]+end.means[,"eph.se"]

ggplot(end.means, aes(x=Sand, y=eph)) + geom_bar(stat="identity", fill="#56B4E9") +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("average surviving Ephedra plant")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") 

## most germinants
maximum.germ <- substrate %>% group_by(Micro, Sand, Rep) %>% summarize(eph.max=max(survival))

max.eph <- maximum.germ %>% group_by(Sand) %>% summarize(eph=mean(eph.max,na.rm=T),eph.se=se(eph.max, na.rm=T)) 
max.eph <- data.frame(max.eph)

max.eph[,"lower"] <- max.eph[,"eph"]-max.eph[,"eph.se"]
max.eph[,"upper"] <- max.eph[,"eph"]+max.eph[,"eph.se"]


ggplot(max.eph, aes(x=Sand, y=eph)) + geom_point(fill="black", size=4) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2)+ylab("total Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("Percentage of sand in soil") + ylim(0,1) + stat_smooth(method="lm", formula= y~poly(x,2),data=max.eph, aes(x=as.numeric(Sand), y=eph), lwd=1, color="black", lty=2)

m1 <- lm(eph~poly(as.numeric(Sand),2), data=max.eph)
summary(m1)

## 
## Call:
## lm(formula = eph ~ poly(as.numeric(Sand), 2), data = max.eph)
## 
## Residuals:
##        1        2        3        4        5 
## -0.01250  0.01875  0.01875 -0.04375  0.01875 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.60625    0.01768  34.295 0.000849 ***
## poly(as.numeric(Sand), 2)1  0.09882    0.03953   2.500 0.129612    
## poly(as.numeric(Sand), 2)2 -0.23385    0.03953  -5.916 0.027402 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03953 on 2 degrees of freedom
## Multiple R-squared:  0.9538, Adjusted R-squared:  0.9075 
## F-statistic: 20.62 on 2 and 2 DF,  p-value: 0.04624

Test of limiting factors - brome

recruit[is.na(recruit)] <- 0

## water
water <- subset(recruit, Treatment == "control" | Treatment=="water")

## germination
water.germ <- water %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(water.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=water, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                          299    1275.05             
## Density      1   849.52       298     425.53   <2e-16 ***
## Lvl          2     4.50       296     421.03   0.1053    
## Density:Lvl  2     0.00       294     421.03   0.9985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## biomass
water.bio <- water %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(water.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## clipping
clip <- subset(recruit, Treatment == "control" | Treatment=="clipped")

## germination
clip.germ <- clip %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(clip.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass
clip.bio <- clip %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(clip.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=clip, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                          299    1332.93             
## Density      1   863.30       298     469.63   <2e-16 ***
## Lvl          2     3.66       296     465.97   0.1605    
## Density:Lvl  2     0.04       294     465.93   0.9814    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## shade
shade <- subset(recruit, Treatment == "control" | Treatment=="shade")

## germination
shade.germ <- shade %>% group_by(Density, Lvl) %>% summarize(eph=mean(ephedra.end ),brome=mean(brome.begin))

ggplot(shade.germ) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("number of Bromus germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass
shade.bio <- shade %>% group_by(Density, Lvl) %>% summarize(eph=mean(Ephedra.biomass ),brome=mean(brome.biomass))

ggplot(shade.bio) + geom_jitter(aes(x=Density, y=brome, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Brome")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=shade, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          299    1221.52              
## Density      1   804.97       298     416.55 < 2.2e-16 ***
## Lvl          2    13.34       296     403.21  0.001267 ** 
## Density:Lvl  2     0.22       294     402.98  0.894507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test of limiting factors - Ephedra

## water
## germination

ggplot(water.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass

ggplot(water.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## clipping

## germination

ggplot(clip.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

## biomass

ggplot(clip.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(ephedra.end~Density *Lvl, data=clip, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                          299     220.71           
## Density      1  2.73641       298     217.97  0.09808 .
## Lvl          2  0.34733       296     217.63  0.84058  
## Density:Lvl  2  0.47792       294     217.15  0.78745  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## shade

## germination

ggplot(shade.germ) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("number of Ephedra germinants")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

##biomass
ggplot(shade.bio) + geom_jitter(aes(x=Density, y=eph, fill=Lvl, color=Lvl), size=2)+ylab("final biomass of Ephedra")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16))

m1 <- glm(brome.begin~Density *Lvl, data=shade, family=poisson)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: brome.begin
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                          299    1221.52              
## Density      1   804.97       298     416.55 < 2.2e-16 ***
## Lvl          2    13.34       296     403.21  0.001267 ** 
## Density:Lvl  2     0.22       294     402.98  0.894507    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Summarized GLMs for manuscript

## Water
m1 <- glm(ephedra.end~ Density *Lvl, data=water, family=binomial)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                          299     397.45           
## Density      1  2.78375       298     394.66  0.09522 .
## Lvl          2  0.20096       296     394.46  0.90440  
## Density:Lvl  2  0.08676       294     394.37  0.95755  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1.bio <- glm(Ephedra.biomass~Density *Lvl, data=subset(water, Ephedra.biomass>0 & Ephedra.biomass<0.6), family=Gamma)
anova(m1.bio, test="Chisq")

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Ephedra.biomass
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                          106     39.539         
## Density      1  0.11078       105     39.428   0.5364
## Lvl          2  0.10032       103     39.328   0.8411
## Density:Lvl  2  0.53107       101     38.796   0.4001

## Shade
m2 <- glm(ephedra.end~ Density *Lvl, data=shade, family=binomial)
anova(m2, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                          299     368.20            
## Density      1   0.7330       298     367.46 0.391927   
## Lvl          2  11.4709       296     355.99 0.003229 **
## Density:Lvl  2   1.4327       294     354.56 0.488533   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2.bio <- glm(Ephedra.biomass~Density *Lvl, data=subset(shade, Ephedra.biomass>0 & Ephedra.biomass < 0.7), family=Gamma)
anova(m2.bio, test="Chisq")

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Ephedra.biomass
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                           67     25.637           
## Density      1  0.39889        66     25.238  0.23231  
## Lvl          2  2.02160        64     23.217  0.02691 *
## Density:Lvl  2  0.02846        62     23.188  0.95038  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Clipping
m3 <- glm(ephedra.end~ Density *Lvl, data=clip, family=binomial)
anova(m3, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                          299     393.19           
## Density      1   4.2054       298     388.98   0.0403 *
## Lvl          2   0.5545       296     388.43   0.7579  
## Density:Lvl  2   0.8190       294     387.61   0.6640  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3.bio <- glm(Ephedra.biomass~Density *Lvl, data=subset(clip, Ephedra.biomass>0 & Ephedra.biomass < 0.7), family=Gamma)
anova(m3.bio, test="Chisq")

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Ephedra.biomass
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                          106     48.613             
## Density      1   0.6174       105     47.996   0.1774    
## Lvl          2   9.2947       103     38.701 1.13e-06 ***
## Density:Lvl  2   0.9905       101     37.711   0.2324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m4 <-  glm(ephedra.end~brome.begin, data=recruit, family=binomial)
anova(m4, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: ephedra.end
## 
## Terms added sequentially (first to last)
## 
## 
##             Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                          699     901.37         
## brome.begin  1   2.0214       698     899.35   0.1551

m4.bio <-  glm(Ephedra.biomass~brome.begin, data=subset(recruit, Ephedra.biomass>0 & Ephedra.biomass < 0.7), family=Gamma)
anova(m4.bio, test="Chisq")

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Ephedra.biomass
## 
## Terms added sequentially (first to last)
## 
## 
##             Df  Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                           209     89.330         
## brome.begin  1 0.0012905       208     89.328    0.953

## brome comparisons
library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:tidyr':
## 
##     expand

## Water
m1 <- glm(brome.begin/Density~ Lvl, weights=Density, data=subset(water, Density>0), family=binomial)
anova(m1, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: brome.begin/Density
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                   239     475.39            
## Lvl   2   10.831       237     464.56 0.004447 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Shade
m2 <- glm(brome.begin/Density~ Lvl, weights=Density, data=subset(shade, Density>0), family=binomial)
anova(m2, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: brome.begin/Density
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   239     429.07              
## Lvl   2   31.496       237     397.57 1.448e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Clipping
m3<- glm(brome.begin/Density~ Lvl, weights=Density, data=subset(clip, Density>0), family=binomial)
anova(m3, test="Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: brome.begin/Density
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                   239     528.32           
## Lvl   2    8.962       237     519.35  0.01132 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Brome density vs Ephedra

dense <- recruit %>% filter(Ephedra.biomass>0) %>%  group_by(brome.begin) %>% summarize(bio=mean(Ephedra.biomass), recruit=mean(ephedra.end), germ=mean(ephedra.begin), abv.grd=mean(Ephedra.above),blw.grd=mean(Ephedra.below))

ggplot(dense, aes(x=brome.begin, y=bio))+geom_point()+ theme_Publication() +ylab("E. californica biomass") + xlab("brome density") + geom_smooth(method=lm, formula= y ~ poly(x,2))

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

m1 <- lm(bio~poly(brome.begin,2), data=dense)

## recruitment no effect
ggplot(dense, aes(x=brome.begin, y=recruit))+geom_point()+ theme_Publication() +ylab("E. californica biomass") + xlab("brome density") #+ geom_smooth(method=lm, formula= y ~ x)

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

m2 <- lm(recruit~brome.begin, data=dense)

## germination no effect
ggplot(dense, aes(x=brome.begin, y=germ))+geom_point()+ theme_Publication() +ylab("E. californica biomass") + xlab("brome density") #+ geom_smooth(method=lm, formula= y ~ x)

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

m3 <- lm(germ~brome.begin, data=dense)


## above ground biomass 
ggplot(dense, aes(x=brome.begin, y=abv.grd))+geom_point()+ theme_Publication() +ylab("above-ground biomass") + xlab("B. madritensis") + geom_smooth(method=lm, formula= y ~ x)

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

m4 <- lm(abv.grd~brome.begin, data=dense)
summary(m4)

## 
## Call:
## lm(formula = abv.grd ~ brome.begin, data = dense)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.010818 -0.003872 -0.001258  0.004408  0.013263 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0365715  0.0034310  10.659 4.21e-08 ***
## brome.begin -0.0009681  0.0003850  -2.515   0.0247 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007263 on 14 degrees of freedom
## Multiple R-squared:  0.3112, Adjusted R-squared:  0.262 
## F-statistic: 6.325 on 1 and 14 DF,  p-value: 0.02474

## below ground biomass
ggplot(dense, aes(x=brome.begin, y=blw.grd))+geom_point()+ theme_Publication() +ylab("below-ground biomass") + xlab("B. madritensis density") + geom_smooth(method=lm, formula= y ~ poly(x,2))

## Warning: `legend.margin` must be specified using `margin()`. For the old
## behavior use legend.spacing

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$label), x$x, x
## $y, : font family not found in Windows font database

m5 <- lm(abv.grd~poly(brome.begin,2), data=dense)
summary(m5)

## 
## Call:
## lm(formula = abv.grd ~ poly(brome.begin, 2), data = dense)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.009539 -0.003013 -0.001133  0.002680  0.011860 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.029250   0.001764  16.586 3.97e-10 ***
## poly(brome.begin, 2)1 -0.018265   0.007054  -2.589   0.0225 *  
## poly(brome.begin, 2)2 -0.009570   0.007054  -1.357   0.1980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.007054 on 13 degrees of freedom
## Multiple R-squared:  0.3966, Adjusted R-squared:  0.3038 
## F-statistic: 4.273 on 2 and 13 DF,  p-value: 0.03748